Hannah L. Owens
28 January, 2020
![]()
![]()


The libraries we're using
library(ape)
#The Swiss army knife of comparative phylogenetics R packages
library(geiger)
#A package for "macroevolutionary simulation and estimating parameters related to diversification"
library(phytools)
#Very useful for visualization particularly, great blog support
istioTree <- read.nexus("Data/Fish_12Tax_time_calibrated.tre")
istioTree
Phylogenetic tree with 12 tips and 11 internal nodes.
Tip labels:
Istiompax_indica, Istiophorus_albicans, Istiophorus_platypterus, Kajikia_albida, Kajikia_audax, Makaira_nigricans, ...
Rooted; includes branch lengths.
str(istioTree)
List of 4
$ edge : int [1:22, 1:2] 13 14 15 16 17 17 16 15 18 19 ...
$ edge.length: num [1:22] 19.98 39.23 6.64 5.47 8.79 ...
$ Nnode : int 11
$ tip.label : chr [1:12] "Istiompax_indica" "Istiophorus_albicans" "Istiophorus_platypterus" "Kajikia_albida" ...
- attr(*, "class")= chr "phylo"
- attr(*, "order")= chr "cladewise"
What does the file look like?
head(istioTree$edge) #Two columns, start node and end node/tip; tips labeled first
[,1] [,2]
[1,] 13 14
[2,] 14 15
[3,] 15 16
[4,] 16 17
[5,] 17 3
[6,] 17 2
head(istioTree$edge.length) #Gives lengths of branches by edge
[1] 19.976147 39.227748 6.635844 5.473333 8.786535 8.786535
What does the file look like?
istioTree$tip.label #Gives a list of tip labels, order corresponds to edge number
[1] "Istiompax_indica" "Istiophorus_albicans"
[3] "Istiophorus_platypterus" "Kajikia_albida"
[5] "Kajikia_audax" "Makaira_nigricans"
[7] "Tetrapturus_angustirostris" "Tetrapturus_belone"
[9] "Tetrapturus_georgii" "Tetrapturus_pfluegeri"
[11] "Trachurus_trachurus" "Xiphias_gladius"
is.ultrametric(istioTree) #Checks to see if the tree is ultrametric
[1] TRUE
The basics
plot(istioTree, cex = 4) #Basic view (phylogram)
nodelabels(cex = 4)
Finding tips that match
colors <- c("black", "red")
matches <- grep("Kajikia", istioTree$tip.label, value = T)
matches
[1] "Kajikia_albida" "Kajikia_audax"
Highlighting matching taxa
plot(istioTree, cex = 4,
tip.col = ifelse(istioTree$tip.label %in% matches,
'red','black'))
Finding the node with a most-recent common ancestor.
#Remember, these are our matches
matches
[1] "Kajikia_albida" "Kajikia_audax"
#What is the node number?
findMRCA(tree = istioTree, tips = matches, type = "node")
[1] 20
We have given you a phylogeny of clupeiform fishes from Bloom and Lovejoy, 2014. Load it into R and display the tree. The file is called “Clupeiform phylogeny.nex”.
Bonus challenge: Which node number represents the most recent common ancestor for the Anchoa clade?
The node of the Anchoa MRCA is 178.
anchovyTree <- read.nexus("Data/Clupeiform phylogeny.nex")
matches <- grep("Anchoa", anchovyTree$tip.label, value = T)
plot(anchovyTree, cex = 0.5,
tip.col = ifelse(anchovyTree$tip.label %in% matches,'red','black'))
findMRCA(tree = anchovyTree, tips = matches, type = "node")
[1] 178
nodelabels(node = 178, text = "MRCA", bg = F, col = "red")
Getting fancy…
plot(istioTree, type = "cladogram", use.edge.length = F,
edge.width = 2, font = 4, cex = 2.5)
Getting fancy…
plot(istioTree, type = "fan", rotate.tree = 30,
main = "Phylogeny of Billfishes", cex.main = 3.5,
label.offset = 2, cex = 2.5,
tip.color = c("blue", "red", "red", "yellow2",
"yellow2", "orange1", "chartreuse4",
"chartreuse4", "chartreuse4",
"chartreuse4", "black", "purple"))
Try to imitate this tree. It is a subset of the anchovy tree trimmed using extract.clade() like this:
subsetAnchovyTree <- extract.clade(anchovyTree, node = 160).
Use ?ape::plot.phylo() for help.
Here's how I made that tree. How'd you do?
subsetAnchovyTree <- extract.clade(anchovyTree, node = 160)
plot(subsetAnchovyTree, type = "fan", rotate.tree = 45, open.angle = 180,
main = "Phylogeny of Anchovies", cex.main = 5, col.main = "green",
label.offset = 2, cex = 2.5, tip.color = c("blue", "red"))
Multiple characters in a .csv file
characterTable <- read.csv("Data/CodingTableThresh95.csv", row.names = 1)
characterTable[,1:4] #Showing the data
bathmin bathmax SSSmin SSS8max
Istiompax_indica -7000 -485 3416 3651
Istiophorus_albicans -6275 -1 3393 3741
Istiophorus_platypterus -7272 -1 3359 3639
Kajikia_albida -5750 -1267 3561 3681
Kajikia_audax -7000 -1 3451 3638
Makaira_nigricans -5750 -1107 3437 3694
Tetrapturus_angustirostris -7000 -841 3425 3649
Tetrapturus_belone -1950 -95 3677 3800
Tetrapturus_georgii -6414 -3392 3458 3741
Tetrapturus_pfluegeri -6693 -2112 3437 3741
Xiphias_gladius -6738 -748 3337 3748
Combining the tree and the character data
treeWData <- treedata(istioTree, characterTable, sort = T)
Warning in treedata(istioTree, characterTable, sort = T): The following tips were not found in 'data' and were dropped from 'phy':
Trachurus_trachurus
#Trims out taxa that are missing from the tree or the character table
Original tree versus trimmed tree
Combining the tree and the character data
treeWData$data[,1:4] #Here's the data
bathmin bathmax SSSmin SSS8max
Istiompax_indica -7000 -485 3416 3651
Istiophorus_albicans -6275 -1 3393 3741
Istiophorus_platypterus -7272 -1 3359 3639
Kajikia_albida -5750 -1267 3561 3681
Kajikia_audax -7000 -1 3451 3638
Makaira_nigricans -5750 -1107 3437 3694
Tetrapturus_angustirostris -7000 -841 3425 3649
Tetrapturus_belone -1950 -95 3677 3800
Tetrapturus_georgii -6414 -3392 3458 3741
Tetrapturus_pfluegeri -6693 -2112 3437 3741
Xiphias_gladius -6738 -748 3337 3748
You've already read in the anchovy phylogeny. Now read in the character table, Bloom_etal_2018_Reduced_Dataset.csv, from Bloom et al., 2018.
How many taxa are represented in the phylogeny and the character dataset?
anchovyCharacters <- read.csv(file = "Data/Bloom_etal_2018_Reduced_Dataset.csv",
row.names = 1)
anchovyTreeWData <- treedata(anchovyTree, anchovyCharacters, sort = T);
plot(anchovyTreeWData$phy, main = "Anchovies!")
paste0("There are ", nrow(anchovyTreeWData$data), " taxa represented in both datasets.")
[1] "There are 50 taxa represented in both datasets."
Plotting discrete traits on tips.
library(RColorBrewer) # Accessory package with better color
#Entering in discrete tip values
discChar <- c(1, 1, 1, 3, 3, 1, 2, 2, 2, 2, 1)
#Plotting
plot(treeWData$phy, main = "Preferred prey", cex.main = 4,
label.offset = 2.5, no.margin = F, cex = 3)
#Setting up colors.
reconCol <- brewer.pal(3, "RdYlBu")
#Plot tips
tiplabels(text = rep(" ", 11), tip = seq(1,11,1), cex = 3,
frame = "circle", bg = reconCol[discChar])
#Plot legend
legend(x = 0, y = 5, fill = reconCol, bty = "n", cex = 3,
legend = c("Sardines", "Pompano", "Hamburgers"))
Plotting discrete traits on tips.
You have a tree united with characters. It's plotting time!
Your challenge: Color-code the tip labels according to whether a species is diadromous or not.
# Set up colors
labelCols <- c("red", "blue");
diadLabs <- labelCols[as.factor(anchovyTreeWData$data[,3])]
#Plotting
plot(anchovyTreeWData$phy, main = "Diadromy", cex.main = 4,
no.margin = F, cex = 3,
label.offset = 2.5)
tiplabels(text = rep(" ", length(anchovyTreeWData$phy$tip.label)),
cex = 1, frame = "circle", bg = diadLabs)
#Plot legend
legend(x = 0, y = 50, bty = "n", cex = 2,
fill = labelCols,
legend = c("Diadromous", "Non-diadromous")
Plotting continuous characters using a color gradient
library(plotrix); #Contains a useful way of producing a color legend
#Set up colors and assigning values
colPal <- brewer.pal(9, "Reds")
normedTips <- round((treeWData$data[,8]-min(treeWData$data[,8]))/
(max(treeWData$data[,8])-min(treeWData$data[,8]))*8,0)
reconCol <- colPal[normedTips+1]
reconCol
[1] "#FB6A4A" "#FC9272" "#FC9272" "#FEE0D2" "#FC9272" "#FEE0D2" "#FB6A4A"
[8] "#FFF5F0" "#FEE0D2" "#FCBBA1" "#67000D"
Plotting continuous characters using a color gradient
# Do the plot
plot(treeWData$phy, main = "Niche breadth", label.offset = 2.5,
no.margin = F)
#Label Tips
tiplabels(text = rep(" ", 11), seq(1,11,1),
cex = 0.65, frame = "circle", bg = reconCol[1:11])
#Making a key for the plot
col.labels<-c("Narrow", "Broad")
color.legend(11,0,13,4, col.labels,
colPal,gradient="y") #Numbers are plot coordinates
Plotting continuous characters using a color gradient
You guessed it. Time for a continuous character plot.
Your challenge: Color-code the tip labels according to log body size.
Bonus challenge: Add another tip label showing trophic position.
(Hint: check out ?tiplabels())
# Set up colors
logBody <- as.numeric(anchovyTreeWData$data[,1])
colPalBody <- rev(brewer.pal(9, "Reds"))
normedTipsB <- round((logBody-min(logBody))/
(max(logBody)-min(logBody))*8,0)
tipColBody <- colPalBody[normedTipsB+1]
trophPos <- as.numeric(anchovyTreeWData$data[,2])
colPalTroph <- rev(brewer.pal(9, "Blues"))
normedTipsT <- round((trophPos-min(trophPos))/
(max(trophPos)-min(trophPos))*8,0)
tipColTroph <- colPalTroph[normedTipsT+1]
#Plotting
plot(anchovyTreeWData$phy, main = "Body Size and Trophic level",
cex.main = 4, no.margin = F, cex = 3, label.offset = 5)
tiplabels(cex = 3, adj = 1, pch = 21, bg = tipColBody)
tiplabels(cex = 3, adj = 3.5, pch = 21, bg = tipColTroph)
#Making keys for the plot
col.labels<- round(c(min(
as.numeric(anchovyTreeWData$data[,1])),
max(as.numeric(anchovyTreeWData$data[,1]))),2)
text(5,46, "Body Size", cex = 3)
color.legend(0,40,4,44,col.labels,rev(colPalBody),
gradient="y", cex = 3, align = "rb")
col.labels<-round(c(min(
as.numeric(anchovyTreeWData$data[,2])),
max(as.numeric(anchovyTreeWData$data[,2]))),2)
text(5,36, "Trophic Level", cex = 3)
color.legend(0,30,4,34,col.labels,rev(colPalTroph),
gradient="y", cex = 3, align = "rb")
library(plotrix) #Accessories for phylogenies
charModel <- fitContinuous(treeWData$phy,
treeWData$data[,8],
model = "BM",
SE = std.error(treeWData$data[,8]));
charModel #How the results look
GEIGER-fitted comparative model of continuous data
fitted 'BM' model parameters:
sigsq = 43607453068390600.000000
z0 = 2554614105.734349
model summary:
log-likelihood = -244.522632
AIC = 493.045264
AICc = 494.545264
free parameters = 2
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
frequency of best fit = 1.00
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
The relevant stats from fitContinuous()
charModel$opt
$sigsq
[1] 4.360745e+16
$z0
[1] 2554614106
$lnL
[1] -244.5226
$method
[1] "Brent"
$k
[1] 2
$aic
[1] 493.0453
$aicc
[1] 494.5453
Run fitContinuous() using two models, Brownian motion and early burst. Which model fits the data better (smallest AICc score).
NOTE: You may run into problems with the $data object if R converted your numbers to characters. Use “as.numeric()” and “names” to fix this!
Oh hey, look at that. Brownian motion is a better fit.
# Data prep
bodySize <- as.numeric(anchovyTreeWData$data[,1])
names(bodySize) <- rownames(anchovyTreeWData$data)
# Fit models
chovyBMcharModel <- fitContinuous(anchovyTreeWData$phy,
bodySize,
model = "BM",
SE = std.error(bodySize))
chovyEBModel <- fitContinuous(anchovyTreeWData$phy,
bodySize,
model = "EB",
SE = std.error(bodySize))
paste0("AICc for BM: ", chovyBMcharModel$opt$aicc, ". AICc for EB: ", chovyEBModel$opt$aicc, ".")
[1] "AICc for BM: 17.1792349761662. AICc for EB: 19.4459022517872."
reconAICcTable <- matrix(data = NA, nrow = 5, ncol = ncol(treeWData$data))
colnames(reconAICcTable) <- colnames(treeWData$data)
rownames(reconAICcTable) <- c("BM", "OU", "EB", "trend", "white")
reconRelWeightTable <- reconAICcTable; #Both tables have the same structure
This loop fits a set of character evolution models for each character.
count <- 1;
while (count <= ncol(treeWData$data)){
reconAICcTable[1,count] <- fitContinuous(treeWData$phy,
treeWData$data[,count], model = "BM",
SE = std.error(treeWData$data[,count]))$opt$aicc
reconAICcTable[2,count] <- fitContinuous(treeWData$phy,
treeWData$data[,count], model = "OU",
SE = std.error(treeWData$data[,count]),
bounds = list(alpha = c(min = exp(-500),
max = exp(10))))$opt$aicc;
reconAICcTable[3,count] <- fitContinuous(treeWData$phy,
treeWData$data[,count], model = "EB",
SE = std.error(treeWData$data[,count]),
bounds = list(a = c(min = -3, max = 3)))$opt$aicc
reconAICcTable[4,count] <- fitContinuous(treeWData$phy,
treeWData$data[,count], model = "trend",
SE = std.error(treeWData$data[,count]))$opt$aicc
reconAICcTable[5,count] <- fitContinuous(treeWData$phy,
treeWData$data[,count], model = "white",
SE = std.error(treeWData$data[,count]))$opt$aicc
deltaAICc <- (reconAICcTable[,count])-min(reconAICcTable[,count])
reconRelWeightTable[,count] <- exp(-.5*deltaAICc)/
sum(exp(-.5*deltaAICc))
count <- count + 1
}
Relative weights of AICc, rounded to four digits:
round(reconRelWeightTable, digits = 5)
bathmin bathmax SSSmin SSS8max SST13min SST13max tempRange nicheRange
BM 0.00094 0.00988 0.00672 0.02128 0.20715 0.00089 0.07440 0.45167
OU 0.10939 0.10807 0.10852 0.12026 0.08281 0.10939 0.10508 0.07090
EB 0.10938 0.10805 0.10852 0.12026 0.08281 0.10939 0.10508 0.07090
trend 0.00038 0.00350 0.00248 0.00709 0.04512 0.00037 0.02010 0.06776
white 0.77991 0.77050 0.77375 0.73112 0.58210 0.77996 0.69534 0.33876
Diet
BM 0.64061
OU 0.11528
EB 0.11528
trend 0.12883
white 0.00000
What is the best model of character evolution for log body size and trophic position in anchovies?
Remember R changed all your numbers to characters. Grr.
Step 1: Set up the data and the table.
# Data
chovyCharTab <- cbind(as.numeric(anchovyTreeWData$data[,1]),as.numeric(anchovyTreeWData$data[,2]))
rownames(chovyCharTab) <- rownames(anchovyTreeWData$data)
colnames(chovyCharTab) <- colnames(anchovyTreeWData$data)[1:2]
# Results tables
reconAICcTable <- matrix(data = NA, nrow = 5, ncol = ncol(chovyCharTab))
colnames(reconAICcTable) <- colnames(chovyCharTab)
rownames(reconAICcTable) <- c("BM", "OU", "EB", "trend", "white")
Step 2: Run that loop!
count <- 1;
while (count <= ncol(chovyCharTab)){
reconAICcTable[1,count] <- fitContinuous(anchovyTreeWData$phy,
chovyCharTab[,count], model = "BM",
SE = std.error(chovyCharTab[,count]))$opt$aicc;
reconAICcTable[2,count] <- fitContinuous(anchovyTreeWData$phy,
chovyCharTab[,count], model = "OU",
SE = std.error(chovyCharTab[,count]),
bounds = list(alpha = c(min = exp(-500),
max = exp(10))))$opt$aicc;
reconAICcTable[3,count] <- fitContinuous(anchovyTreeWData$phy,
chovyCharTab[,count], model = "EB",
SE = std.error(chovyCharTab[,count]),
bounds = list(a = c(min = -3, max = 3)))$opt$aicc
reconAICcTable[4,count] <- fitContinuous(anchovyTreeWData$phy,
chovyCharTab[,count], model = "trend",
SE = std.error(chovyCharTab[,count]))$opt$aicc
reconAICcTable[5,count] <- fitContinuous(anchovyTreeWData$phy,
chovyCharTab[,count], model = "white",
SE = std.error(chovyCharTab[,count]))$opt$aicc
count <- count + 1
}
Step 3: Which model?
reconAICcTable
logbodysize trophic_position
BM 17.17923 -100.7962
OU 17.31522 -156.8098
EB 17.31522 -156.6843
trend 16.98079 -107.3646
white 49.82883 -159.0762
phytools::phylosig(treeWData$phy, x=treeWData$data[,1],method = "K")
[1] 0.1342713
sigTable <- matrix(data = NA, ncol = 4, nrow = ncol(treeWData$data))
rownames(sigTable) <- colnames(treeWData$data)
colnames(sigTable) <- c("Lambda", "Lambda P-value", "K", "K P-value")
Pagel's Lambda, Blomberg's K, and their respective P-values for each character
count <- 1
while(count <= ncol(treeWData$data)){
temp <- phytools::phylosig(treeWData$phy, treeWData$data[,count],
method = "lambda", test = T); #Calculate Pagel's lambda
sigTable[count,1] <- temp$lambda
sigTable[count,2] <- temp$P
temp <- phytools::phylosig(treeWData$phy, treeWData$data[,count],
method = "K", test = T); #Calculate Blomberg's k
sigTable[count,3] <- temp$K
sigTable[count,4] <- temp$P
count <- count + 1
}
Et voila!
sigTable
Lambda Lambda P-value K K P-value
bathmin 6.901244e-05 1.000000e+00 0.1342713 0.835
bathmax 6.901244e-05 1.000000e+00 0.1958716 0.700
SSSmin 6.754471e-05 1.000000e+00 0.2386497 0.558
SSS8max 6.901244e-05 1.000000e+00 0.2418311 0.546
SST13min 6.574920e-01 5.737320e-01 0.4547630 0.100
SST13max 6.901244e-05 1.000000e+00 0.1375421 0.820
tempRange 5.217411e-05 1.000000e+00 0.4011244 0.163
nicheRange 7.364686e-01 1.068847e-01 0.8079658 0.060
Diet 1.043882e+00 7.949707e-05 1.1791841 0.002
What is the phylogenetic signal for log body size and trophic position in anchovies?
Step 1: Set up table for results.
#Set up a table for results
sigTable <- matrix(data = NA, ncol = 4, nrow = ncol(chovyCharTab))
rownames(sigTable) <- colnames(chovyCharTab)
colnames(sigTable) <- c("Lambda", "Lambda P-value", "K", "K P-value")
Step 2: Loop through each column.
#Loop to calculate Blomberg's K and P-value for each character
count <- 1
while(count <= ncol(chovyCharTab)){
temp <- phytools::phylosig(anchovyTreeWData$phy, chovyCharTab[,count],
method = "lambda", test = T); #Calculate Pagel's lambda
sigTable[count,1] <- temp$lambda
sigTable[count,2] <- temp$P
temp <- phytools::phylosig(anchovyTreeWData$phy, chovyCharTab[,count],
method = "K", test = T) #Calculate Blomberg's k
sigTable[count,3] <- temp$K
sigTable[count,4] <- temp$P
count <- count + 1
}
Step 3: How'd we do?
reconAICcTable
logbodysize trophic_position
BM 17.17923 -100.7962
OU 17.31522 -156.8098
EB 17.31522 -156.6843
trend 16.98079 -107.3646
white 49.82883 -159.0762
sigTable
Lambda Lambda P-value K K P-value
logbodysize 9.954688e-01 9.920151e-09 0.7099999 0.001
trophic_position 5.186526e-05 1.000000e+00 0.1009663 0.809
recon <- ace(x = treeWData$data[,8], phy = treeWData$phy,
type = "continuous", method = "GLS",
corStruct = corBrownian(value=1, treeWData$phy))
recon
Ancestral Character Estimation
Call: ace(x = treeWData$data[, 8], phy = treeWData$phy, type = "continuous",
method = "GLS", corStruct = corBrownian(value = 1, treeWData$phy))
$ace
12 13 14 15 16 17 18
2589649968 1428149517 1272412103 1405633461 1393716796 1296733993 1235410112
19 20 21
1084046522 1063892522 721486716
$CI95
[,1] [,2]
12 2589649968 2589649968
13 1428149513 1428149522
14 1272412099 1272412107
15 1405633458 1405633465
16 1393716792 1393716799
17 1296733990 1296733997
18 1235410110 1235410114
19 1084046519 1084046525
20 1063892520 1063892525
21 721486714 721486719
The setup
reconMatrix <- matrix(nrow=(length(treeWData$phy$tip.label) + treeWData$phy$Nnode),
ncol=1)
colnames(reconMatrix) <- colnames(treeWData$data)[18]
rownames(reconMatrix) <- c(treeWData$phy$tip.label,
seq(from = length(unlist(treeWData$phy$tip.label)) + 1,
to = nrow(reconMatrix)))
Reconstruct niche breadth
recon <- ace(treeWData$data[,8], treeWData$phy, type = "continuous",
method= "GLS", corStruct = corBrownian(value=1, treeWData$phy))
reconMatrix[,1] <- round(c(treeWData$data[,8], recon$ace), 2)
#Putting data from reconstructions into a matrix.
reconMatrix <-as.matrix(reconMatrix)
Plotting reconstructed niche breadth
library(RColorBrewer) # Accessory package with better color
colPal <- brewer.pal(9, "RdYlBu")
reconCol <- colPal[round((reconMatrix[,1]-min(reconMatrix[,1]))/(max(reconMatrix[,1])-min(reconMatrix[,1]))*8,0)+1]; #Determining what color will be displayed at each node and tip.
plot(treeWData$phy, main = "Reconstruction of niche breadth", label.offset = 2.5, no.margin = F) #Set down the plot.
nodelabels(text = rep(" ", 21-12+1), node = seq(12,21,1), cex = 0.65, frame = "circle", bg = reconCol[12:21]) #Label the nodes.
tiplabels(text = rep(" ", 11), seq(1,11,1), cex = 0.65, frame = "circle", bg = reconCol[1:11]) #Label the tips.
#Making a key for the plot
col.labels<-c("Narrow", "Broad")
color.legend(11,0,13,4,col.labels,colPal,gradient="y")
Plotting reconstructed niche breadth
Reconstruct and plot ancestral anchovy body size. Remember that the Ornstein Uhlenbeck evolutionary model was the best fit for this data. So the correlation structure you want to use is corMartins instead of corBrownian.
Step 1: Do the reconstruction.
recon <- ace(chovyCharTab[,1], anchovyTreeWData$phy,
type = "continuous",
method= "GLS", corStruct = corBrownian(value=1, anchovyTreeWData$phy))
reconMatrix <- round(c(chovyCharTab[,1], recon$ace), 2)
#Putting data from reconstructions into a matrix.
reconMatrix <-as.matrix(reconMatrix)
Step 2: Plot it.
colPal <- brewer.pal(9, "PuOr")
reconCol <- colPal[round((reconMatrix[,1]-min(reconMatrix[,1]))/(max(reconMatrix[,1])-min(reconMatrix[,1]))*8,0)+1] #Determining what color will be displayed at each node and tip.
plot(anchovyTreeWData$phy, main = "Reconstruction of body size", label.offset = 2.5, no.margin = F) #Set down the plot.
nodelabels(text = rep(" ", 49), node = seq(51,99,1), cex = 0.65, frame = "circle", bg = reconCol[51:99]) #Label the nodes.
tiplabels(text = rep(" ", 50), seq(1,50,1), cex = 0.65, frame = "circle", bg = reconCol[1:50]) #Label the tips.
#Making a key for the plot
col.labels<-c("Small", "Large")
color.legend(0,50,10,40,col.labels,colPal,gradient="y")
Step 2: Plot it.
The basic function
discChar <- c(1, 1, 1, 3, 3, 1, 2, 2, 2, 2, 1)
#Entering in discrete tip values
recon <- ace(discChar, treeWData$phy, type = "discrete", method = "ML", model = "ER");
#Do reconstruction (equal likelihood for all transitions)
recon #What it looks like
Ancestral Character Estimation
Call: ace(x = discChar, phy = treeWData$phy, type = "discrete", method = "ML",
model = "ER")
Log-likelihood: -7.502645
Rate index matrix:
1 2 3
1 . 1 1
2 1 . 1
3 1 1 .
Parameter estimates:
rate index estimate std-err
1 0.012 0.0081
Scaled likelihoods at the root (type '...$lik.anc' to get them for all nodes):
1 2 3
0.5242880 0.2443479 0.2313642
What the ancestral state reconstructions look like.
recon$lik.anc
1 2 3
[1,] 0.5242879674 0.244347855 2.313642e-01
[2,] 0.7988762539 0.129944244 7.117950e-02
[3,] 0.9653742944 0.021355494 1.327021e-02
[4,] 0.9927999747 0.004170209 3.029816e-03
[5,] 0.6691495041 0.223404708 1.074458e-01
[6,] 0.4347414214 0.386702471 1.785561e-01
[7,] 0.0033931287 0.003106903 9.935000e-01
[8,] 0.0292836926 0.957222775 1.349353e-02
[9,] 0.0014809269 0.997760499 7.585745e-04
[10,] 0.0001420902 0.999761557 9.635278e-05
Plotting the reconstruction results.
plot(treeWData$phy, main = "Reconstruction of preferred prey", cex.main = 5, label.offset = 2.5, no.margin = F, cex = 4)
#Set down the plot
reconCol <- brewer.pal(3, "RdYlBu")
#Setting up colors.
nodelabels(text = rep(" ", 21-12+1), node = seq(12,21,1), pie = recon$lik.anc,
piecol = c(reconCol[1], reconCol[2], reconCol[3]), cex = .9)
#Label the nodes.
tiplabels(text = rep(" ", 11), seq(1,11,1), cex = 3, frame = "circle", bg = reconCol[discChar])
#Label the tips.
legend(x = 0, y = 5, legend = c("Sardines", "Pompano", "Hamburgers"), fill = reconCol, bty = "n", cex = 5)
#Plot the legend.
Plotting the reconstruction results.
Was the most recent common ancestor of anchovies more likely diadromous or not based on our data?
Hint: Use as.factor and as.numeric to transform your tip data from character strings to 1s and 2s!
Step 1: Do the reconstruction.
#Set up the data
diadromy <- as.numeric(as.factor(anchovyTreeWData$data[,"Reg"]))
names(diadromy) <- rownames(anchovyTreeWData$data)
recon <- ace(diadromy, anchovyTreeWData$phy, type = "discrete", method = "ML", model = "ER")
#Do reconstruction (equal likelihood for all transitions)
recon #What it looks like
Ancestral Character Estimation
Call: ace(x = diadromy, phy = anchovyTreeWData$phy, type = "discrete",
method = "ML", model = "ER")
Log-likelihood: -23.89524
Rate index matrix:
1 2
1 . 1
2 1 .
Parameter estimates:
rate index estimate std-err
1 0.0056 0.0018
Scaled likelihoods at the root (type '...$lik.anc' to get them for all nodes):
1 2
0.0991219 0.9008781
Step 2: Plot it.
Plotting the reconstruction results.
plot(anchovyTreeWData$phy, main = "Evolution of Diadromy", cex.main = 5, label.offset = 2.5, no.margin = F, cex = 4)
#Set down the plot
reconCol <- c("Red", "Blue")
#Setting up colors.
nodelabels(text = rep(" ", 49), node = seq(51,99,1),
pie = recon$lik.anc,
piecol = c(reconCol[1], reconCol[2]), cex = .9)
#Label the nodes.
tiplabels(text = rep(" ", 50), seq(1,50,1),
cex = 3, frame = "circle", bg = reconCol[diadromy])
#Label the tips.
legend(x = -6, y = 50, legend = c("Diadromous", "Not"),
fill = reconCol, bty = "n", cex = 3)
#Plot the legend.
Step 2: Plot it.
Plotting the reconstruction results.
Analysis of Phylogenetics and Evolution with R
Phytools